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We propose a model of a nonlinear double- well potential (NDWP), alias a double-well pseudopo- 
tential, with the objective to study an alternative implementation of the spontaneous symmetry 
breaking (SSB) in Bose-Einstein condensates (BECs) and optical media, under the action of a po- 
tential with two symmetric minima. In the limit case when the NDWP structure is induced by the 
local nonlinearity coefficient represented by a set of two delta-functions, a fully analytical solution 
is obtained for symmetric, antisymmetric and asymmetric states. In this solvable model, the SSB 
bifurcation has a fully subcritical character. Numerical analysis, based on both direct simulations 
and computation of stability eigenvalues, demonstrates that, while the symmetric states are sta- 
ble up to the SSB bifurcation point, both symmetric and emerging asymmetric states, as well as 
all antisymmetric ones, are unstable in the model with the delta-functions. In the general model 
with a finite width of the nonlinear-potential wells, the asymmetric states quickly become stable, 
simultaneously with the switch of the SSB bifurcation from the subcritical to supercritical type. 
Antisymmetric solutions may also get stabilized in the NDWP structure of the general type, which 
gives rise to a bistability between them and asymmetric states. The symmetric states require a finite 
norm for their existence, an explanation to which is given. A full diagram for the existence and 
stability of the trapped states in the model is produced. Experimental observation of the predicted 
effects should be possible in BEC formed by several hundred atoms. 

PACS numbers: 03.75.Lm, 03.75.Kk, 05.45.Yv, 42.65.Tg 



I. INTRODUCTION 



The one-dimensional (ID) Schrodinger equation including a symmetric potential structure produces single-particle 
wave functions of a definite parity, even or odd, with the ground state always corresponding to an even function 
without zeros. However, a spatially symmetric Hamiltonian of an interacting many-particle system can give rise to 
asymmetric states, which may be considered as a spontaneous-symmetry-breaking (SSB) effect. At the classical level, 
the SSB occurs in optics, as a result of the interplay between the nonlinearity and waveguiding structures, when the 
strong nonlinearity partly suppresses the linear coupling between parallel guiding cores. In particular, it was shown 
that a stable trapped mode may be asymmetric in a channel waveguide embedded in the self-focusing Kerr medium 
[l[ . The onset of a sharp symmetry-breaking instability in a double- hump two-component spatial optical soliton was 
demonstrated experimentally in a planar nonlinear waveguide . 

A natural setting in which SSB phenomena may arise in the context of nonlinear optics and Bose-Einstein conden- 
sation (BEC) is provided by double- well potentials (DWPs). In the experiment, an effective optical DWP was created 
by a specially designed illumination pattern applied, in the ordinary polarization, to a photorefractive crystal (the 
SSB was observed in a beam with extraordinary polarization, shone through this structure) It was also proposed 
to realize similar effective potentials in coupled nonlinear microcavities and in a structured core of an optical 
fiber ^ . A specific variety of the optical SSB was studied in a model of two parallel-coupled antiwaveguides with the 
self- focusing nonlinearity, which corresponds to an effective double-barrier potential, rather than DWP [Q]. 

Well-known dual-core optical fibers Q , which may serve as a basis for the power-controlled all-optical switching, if 
the Kerr nonlinearity is taken into regard [sj, may also be considered as DWP structures, with the difference that the 
tunneling between two potential wells is replaced by the linear coupling between the cores. In addition to the SSB of 
continuous- wave states [9| , the formation of asymmetric solitons in dual-core fibers was studied in detail theoretically 
[lo| . Similar analysis of the SSB for soliton modes was performed in models of dual-core fiber Bragg gratings with 
the Kerr nonlinearity pd| . and coupled waveguides with the quadratic fl3| and cubic-quintic (Tsj nonlinear terms, 
including a system of linearly coupled complex Ginzburg-Landau equations of the cubic-quintic type (l4| . In Refs. 
[m , [iBl and [l3| , the analysis of the SSB was extended to three-core linearly coupled triangular configurations - for 
optical fibers, Bragg gratings, and complex Ginzburg-Landau equations, respectively. 

The concept of SSB also plays an important role in understanding experimental phenomena in BEC, because, if 
interactions between atoms are strong enough, the ground state of the condensate may not follow the symmetry of 
the trapping potential [18] . In particular, manifestations of SSB were observed in a quenched ferromagnetic state 



2 



of a spinor (three-component) condensate [l9|. In the single-component BEC, a natural setting for the realization 
of SSB may again be provided by DWP configurations. An effectively one-dimensional DWP structure was realized 
experimentally in Ref. [13] ■ Loading a condensate of ^^Rb atoms with the repulsive interaction between them into 
this structure made it possible to observe Josephson oscillations for a small number of atoms, and the macroscopic 
quantum self-trapping featuring an imbalance between populations of the two wells, for a larger number. Parallel to 
the experimental work, numerous theoretical studies of matter-wave DWP settings have been performed, for the cases 
of both repulsion and attraction between atoms. These studies addressed problems such as finite-mode reductions 
[2H (including two-component mixtures [23 )j obtaining analytical results for specific shapes of the potential [1^, 
quantum effects , and some others. Recently investigated tunneling between vortex and antivortex states in EEC 
trapped in a 2D anisotropic potential [1^ belongs to this category too. 

Theoretical analysis was also performed for 2D and 3D extensions of the DWP settings in BEC, which add one 
or two extra dimensions to the model, either without an additional potential, or with a periodic optical-lattice (OL) 
potential acting in these directions. These settings may be approximated, similar to the above-mentioned standard 
model of dual-core optical fibers, by a system of linearly coupled ID [2I] or 2D [23] equations. In a more accurate form, 
nearly-lD solitons can be found as solutions to the full 2D equation that includes the DWP (the potential depends on 
the transverse coordinate, x, allowing solitons to self-trap in the free longitudinal direction, y) [28]. The latter model 
is relevant to the case of the self- attractive nonlinearity. In the case of self-repulsion, dual-core gap solitons have been 
be predicted in the setting with the OL potential applied along direction y [29|. Note that, in any setting, gap solitons 
cannot realize the ground state of the respective system, but, nevertheless, they represent stable configurations, that 
have been created in the experiment using the condensate of ®^Rb atoms with the repulsion between them [s^. 

A general principle, upheld by the analysis in various settings (in nonlinear optics and BEC alike), is that the 
SSB occurs through bifurcations of symmetric or antisymmetric states, in the models with the self-attraction and 
self- repulsion, respectively. As mentioned above, models of the DWP/double-core type, combining cubic attractive 
and quintic repulsive nonlinearities, were studied too [H, [ij, lU [s^]. In the latter case, the competition between 
the self-focusing and self-defocusing against the backdrop of the DWP structure gives rise to specific SSB bifurcation 
diagrams, in the form of non-convex closed loops 's^l, as well as to specific dynamical switching regimes (3l|. Also 
predicted were manifestations of the SSB in a two-component BEC mixture trapped in the DWP structure |22l |. for 
both cases of the self-attraction and self-repulsion. 

All the extensive work on the SSB outlined above was performed in settings based on usual linear potentials of 
the double-well type. The objective of the present work is to propose another physical framework, in which the SSB 
can be predicted in an effective nonlinear double- well potential (NDWP), induced through a spatial modulation of 
the local nonlinearity coefficient. Following the terminology commonly adopted in the solid-state theory [ssj, this 
nonlinear ingredient of the physical model may also be a called a pseudopotential. 

In BEC settings, a pseudopotential structure may be readily induced through spatial modulation of the local value 
of the s-wave scattering length, as{x), which determines the effective BEC nonlinearity. The modulation can be 
implemented, through the Feshbach resonance, by means of a spatially inhomogeneous dc inagnetic field [13], or by a 
resonant optical field, as predicted in Ref. [3 51 ] and demonstrated experimentally in Ref. [36"]. It was also proposed 
to control the Feshbach resonance by dint of dc electric field |33], which can be easily made inhomogeneous too. The 
attractive and repulsive interactions between atoms correspond to < and > 0, respectively; both signs, as well 
as sign-changing patterns, can be used to engineer effective nonlinear potentials. 

So designed pseudopotential lattices have attracted much interest in studies of BEC. In the ID geometry, solitons, 
extended wave patterns, and various dynamical states supported by such structures were studied theoretically [38l.l39| 
(a random nonlinear lattice [40| and pseudopotentials generated by a spatially monotonous ramp of the local scattering 
length [4l| were explored too). Recently, similar states were also considered in nonlinear optics, assuming a periodic 
modulation of the local Kerr coefficient 42] . Some (but much fewer) results were obtained too for 2D settings (4^ . 

However, to the best of our knowledge, SSB phenomena in nonlinear pseudopotentials have not been studied yet. 
In this work, we focus on such effects in NDWP settings, which can be engineered by means of techniques mentioned 
above, using attractive interactions between atoms in BEC, or the self-focusing nonlinearity in optics, as briefiy 
described below. In Section II, we formulate the model and give estimates of characteristic values of related physical 
parameters. Section III reports full analytical solutions corresponding to symmetric, antisymmetric, and asymmetric 
states trapped by the NDWP, in the limit case with the modulation of the local nonlinearity coefficient is represented 
by a set of two Dirac's delta- functions. The relevance of the latter model is stressed, in particular, by the recently 
introduced [39j BEC model with a periodic nonlinear potential of the Kronig-Penney type, whose simplest version 
reduces to a periodic array of delta- functions (SSB effects were not studied in Ref. [33]). 

In section IV, we present numerical results for the general model, in which the delta-functions are replaced by a 
pair of Gaussians of a finite width. In that case, the trapped states are found in a numerical form, and their stability 
is studied by means of direct simulations of slightly perturbed stationary states and also, independently, through 
computation of respective stability eigenvalues for small perturbations. The result is that asymmetric states, which 
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are unstable in the delta- function limit, can be readily stabilized in the general model. In addition, antisymmetric 
states may be stabilized too, in two disjoint regions of the parameter state, giving rise to a bistability involving 
antisymmetric and asymmetric states. A condition for the existence of symmetric states is that their norm must 
exceed a certain threshold value, in terms of their norm, an explanation to which is given. The existence and stability 
of all states, including a line of the SSB bifurcation, are summarized in a single diagram, which is presented in Section 
IV too. Results reported in this paper and perspectives for further work are summarized in Section V. 

II. THE MODEL 

The underlying 3D Gross-Pitaevskii equation for the mean- field wave function, '^{X, Y, Z, T), is taken in the ordinary 
form corresponding to the nearly-lD trap (43 |: 



(1) 



2m 2 m 

with TO the atomic mass and wj^ the frequency providing for the tight confinement of the condensate in the direction 
of i? = \/Y'^ -f Z^. As said above, we assume an axial modulation (and negative sign) of the scattering length, in the 
form of a pair of Gaussians, each of width I and amplitude Aq, with centers set at points X — ±A: 



as{X) 



exp 



(X + Af 

^2 



exp 



(2) 



In what follows below, we measure the axial coordinate in units of A, and, accordingly, time in units of mh? /h, i.e., 
we define 



X EE X/l, a = l/A,t = hT/ (toA^) 



(3) 



Further, following the usual approach to the derivation of the effective ID equation, we take 3D field as a product of 
a slowly varying axial function, ijj{X,T), and the ground-state wave function in the transverse plane, 



* = 



1 



A/27r5/2a^oA2 



i?2 

^1 



exp ( iuj±_T - ^ ) ip{x, t), 



(4) 



with a\ = hl (rabj±) [the scaling factor here is chosen so as to maintain normalization condition ([7]), see below]. 
The substitution of expression ^ in Eq. ^ and averaging in the transverse plane lead to the following ID equation: 



i->Pt = --ipxx + g{x)\i)\^ip, 



(5) 



9{x) 



exp 



(x+ir 



exp 



{x~lf 



where the modulated nonlinearity coefficient is subject to the following normalization condition: 

g{x)dx = 2. 



(6) 



(7) 



Profiles of modulation function which keeps the double- well structure for a < \/2, are shown, for different values 
of a, in Fig. [TJ 

Equation ([5|) conserves two dynamical invariants, viz., the norm and energy (Hamiltonian), 



N 



1 f°° / \ 



dx. 



(8) 



As follows from Eqs. ^ and ([3]), the total number of atoms in the condensate, Af, is related to ID norm N as follows: 

„2 



l"^ {X,Y,Z)\^ dXdYdZ = 



27r3/2a^oA 



N. 



(9) 
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FIG. 1: (Color online) Shapes of double-well function ((6)l, normalized to its maximum value, are shown for different values of 
scaled width a of the individual well. 



For physical parameters relevant to experiments with the condensate of ^Li atoms pa j. i.e., a± ^ 2 /im, ~ 0.5 nm, 
and A ~ 20 /im, characteristic values of the number of atoms in various patterns reported below [see, in particular. Figs. 
[Sljd) and [9] fall into the range of M between ^ 200 and 1000, which is quite sufficient for experimental manipulations 
and observation of the patterns. In the same range of physical parameters, t = 1 in Eq. ([5]), is estimated as being 
tantamount to ^ 10 ms, hence typical time scales for the instability development or intrinsic oscillations of breathers 
induced by the instabilities, that are reported below, are expected to be in the range of 0.1 to 1 s, which is realistic 
to the currently available experimental techniques [4a |. 

In terms of optical settings, a set of two narrow (of width I ~ 1 /^m) parallel stripes with strong local nonlinearity 
can be built, in a planar waveguide, by means of available nanotechnological methods. In that case, the power of the 
laser beam necessary for the self-trapping of transverse nonlinear patterns in the waveguide made of silica may be 
~ 500 kW, \id\ while, using AlGaAs, one may reduce the necessary power to the level of 1 kW 47|. In these settings, 
the characteristic evolution length of the spatial beam can be made shorter than 1 mm. Obviously, transitions between 
states of different types reported in this paper may be relevant to the design of power-controlled optical-switching 
schemes. On the other hand, the description of the planar waveguide with the pair of embedded stripes may require 
a model more general than the one studied here, as it will plausibly combine the transverse modulation of the local 
nonlinearity with a similar linear potential (which is briefly described at the end of the next section) , as the material 
difference between the stripes and host medium ought to affect the linear index of refraction too. 

Stationary localized solutions to Eq. ([5]) are sought for as ip = e^*^*(/)(x), where the chemical potential is negative, 
/i < 0, and real fimction (t>{x) satisfies equation 



2 avTT 



exp 



exp 



= 0. 



(10) 



An analytically solvable version of the model corresponds to the limit of a ^ 0, with Eq. ([5]) going over into 

ii,t = -(1/2)^.. - \b{x + 1) + b(x - 1)] IV'PV', (11) 

where 5(x) is the Dirac's delta-function. Note that rescaled equation pd|) contains no free parameters; however, the 
norm of the solution will play the role of an intrinsic parameter, see below. In the same limit, stationary equation 
takes the form of 



li<l> + (1/2)0" + [(5(.T + 1) + 8{x - 1)] <i? = 0. 



(12) 
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III. THE MODEL WITH THE DELTA-FUNCTIONS: ANALYTICAL SOLUTIONS 
A. Symmetric, antisymmetric, and asymmetric states 
Off points X = ±1, Eq. (fT2|) is linear. A general solution to this equation, decaying at |a;| oo, can be written as 



BieV^\i'^+^) , at a; < -1, 
(x) = { Aoe-%/2H(^-i) + BoeV2H(^+i)^ at - 1< a: < +1, 
^ie-V2R(^-i), ^^^^ 



(13) 



with constant amplitudes Aq, Ai and Bq, Bi. The continuity of the wave function at a; = ±1 imposes two relations on 
them, Bi = Bo + Aoe'^V^^^, Ai = Aq + Bqc^V^^^ , which allows one to eliminate Aq and Bq in favor of Ai and Bi, 



eW^MBi-Ai ^ eV2|MUi-Bi 
Aq — 1== , Bq — . 



sV2|Ml „ 1 



(14) 



Further, the integration of Eq. ()12|) in infinitesimal vicinities of points x — ±1 yields expressions for jumps (A) of 
the first derivative at these points, A (0') |a;=±i — —2 {<p\x=±i) ■ The substitution of solution in these relations 
leads to a system of cubic equations for the amplitudes. 



/H72(i3i-So + AoeV2|Mlj = Bl 
/H72 (^1 - Ao + Boe'V^) = Al 



(15) 
(16) 

After the substitution of expressions (fT4|) into Eqs. (fTH]) and (fTS]) . we end up with two coupled cubic equations for Ai 
and Bi : 



=4V2|mI 



-IB 



(17) 



/^e^VJH (^^/^\A^ - Bi) = (e^V^ - l) A\. 



Solving Eqs. (I17p and (|18p . we first find symmetric and antisymmetric solutions. 



= Bi = A, 



±1 



l + e 



-2V2|mI 



(18) 



(19) 



^1 — —B\ — Aantisym — ± A 



The norm of these solutions, defined as per Eq. ([8]), is 



/2H 



1 _ p-2V2|Ml 



(20) 



1 - e-V2|Ml ± 4^/2Re- V2|mI 



sym.antisym 



(21) 



with + and — corresponding to the symmetric and antisymmetric states, respectively. In the limit of /i — > — oo, 
both expressions (fT9|) and ((20|) yield ^sym (^ = — cx)) = ^antisym (a^ ~ — oo) = ^2|^|. In this Hmit, norm ([2T|) of 
both solutions takes a common value, iVsym,antisym (m — — oo) = 2. On the other hand, in the limit of /i ^ —0, the 
amplitude of the symmetric state vanishes, y4.sym(/i — > —0) « (|/i|/2)^^*, and its norm takes a finite limit value, 
NgymilJ' = 0) = 1/2, while the amplitude of the antisymmetric state remains finite, Aantisyin(M ~^ ^0) = Vv^, and 
its norm diverges, iVantisym(M ^ ^0) « 1/ ^2a/2|/i|^ . It can be easily checked that the decrease of A'gym and increase 
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of A^antisym with the variation of /x from —oo to are not monotonous: iVsyni attains a maximum, (-^sym),jjax ~ 2.08, 
at /X « -1.40, and A^antisym has a minimum, (A^^antisym),!,;!! ~ l-^^' at « -0.58. 

Actually, the finite minimum norm (threshold), necessary for the existence of the symmetric states, which is 
NsymifJ- = 0) = (-^sym)j„in — 1/2 in thc prcsent case, is a generic property, shared by the model with finite width a 
of the nonlinear-potential wells, as shown in detail below in Fig. [9l On the other hand, the existence of the above- 
mentioned maximum value of the norm for the symmetric states, {Nsym)^^,^ ~ 2.08, is a specific feature of the model 
with a = (the one based on the delta- functions) : it is easy to see that, in the model with finite a, the norm of any 
state grows ~ at /i — > —oo. However, this specific feature, which, as a matter of fact, indicates degeneracy of 

the model with a = 0, is less significant for the physical applications, as all symmetric states, for a = and finite a 
alike, are unstable when the norm exceeds its value at thc SSB bifurcation point [further details can be seen in Eq. 
([23| and Fig. [9] below]. 

A point of the SSB bifurcation, which gives rise to a pair of asymmetric solutions splitting off from the symmetric 
one, can be easily found. Indeed, thc symmetry breaking means that symmetric solution (I19|) acquires an infinitesimal 
antisymmetric addition, SAi — —5Bi = 5A. Thus, infinitely close to the bifurcation point, the relevant solution is 
sought for as Ai = Agym + 5 A, Bi = Agym — SA. The substitution of this in Eqs. p?|l and (fT8|) and linearization in 
infinitesimal SA lead to a simple equation that predicts the value of the chemical potential at the bifurcation point: 

exp ^\/2 l/ibif l) = V2, or 

/ibif = - (In 2f /8 « -0.06 . (22) 



At this point, thc amplitude of symmetric solution (fT9|) is Aut — \J (In 2) /3 « 0.481, and the value of norm ([8]), with 
the upper sign, is 

iVbif = 2/3+ (8/27) (3/4 -I- In 2) w 1.09. (23) 

Actually, value (|22p of the chemical potential at the bifurcation point ins approximately the same in thc model with 
finite- width nonlinear-potential wells, up to a ~ 1, as seen from Fig. [3]Jd) presented below. 

The same analysis shows that antisymmetric solution (|20p never gives rise to an antisymmetry-breaking bifurcation. 
Indeed, for this solution the antisymmetry would be broken by an infinitesimal symmetric variation, 6A\ = bB\ = &A, 
i.e., infinitely close to the bifurcation point, the solution would be A\ = Aantisym + ^A, B\ = — ^antisym + ^A. 
Subsequent substitution in Eqs. p7|) and (jl8p and the linearization in 8 A yield an equation for \i that has no real 
solutions. 

Equations p?]) and can be solved analytically for asymmetric states too: 



{^1. 5l}asy.„ - ^ ; ^ • (24) 



Note that full solution (I24p predicts exactly the same bifurcation point as Eq. (|22p . i.e., exp |^-\/2|/ibif | j = \/2 [at this 

point, the second radical in Eq. (|24p vanishes]. These solutions are characterized by the asymmetry ratio, which is 
defined as 

Typical examples of symmetric, asymmetric and antisymmetric states produced by he above analytical solutions are 
displayed in Fig. [2l 



B. Bifurcation diagrams and stability 



The analytical solution given by Eqs. p3|) . p4)) and ((24)) make it possible to plot the bifurcation diagrams in the 
planes of (/x, 8) and {N, 0), which are represented by curves pertaining to a = in Figs. [3]^a,b,c). To generate the 
diagrams, partial norms N± in expression ([25]) for the asymmetric solutions were computed numerically [analytical 
expressions for them are available, but they are very messy, cf. Eqs. (I2ip for the symmetric and antisymmetric states]. 
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FIG. 2: (Color online) Profiles of symmetric, asymmetric, and antisymmetric stationary states in the analytically solvable model 
with the delta-functions, all pertaining to ^ = —0.1. The respective values of the norm are A^symm ~ 1.27, A^asymm ~ 1.07, and 
A'^antisymm ~ 2.18. The mutually symmetric lines (blue and red ones, in the color version of the figure) represent two asymmetric 
states that are mirror images of each other. 

A salient peculiarity of the SSB bifurcation for a = 0, evident in Fig. [Sjjc), is its subcritical character, which means 
that the branches of asymmetric solutions emerge at the bifurcation point as unstable ones, and go in the backward 
direction. A subcritical bifurcation also occurs in the above-mentioned model of the dual-core nonlinear fibers p^ . 
but in that case the asymmetric branches quickly turn in the forward direction, getting stabilized at the turning point. 
A remarkable feature of the present model with a = is that this does not happen, i.e., the bifurcation in this model 
may be called a "fully backward" one: the branches of the asymmetric solutions keep going backward up to the limit 
of O = 1, which corresponds to the asymmetric solutions with /i = — oo (and N = 1, as shown in the following 
subsection). Indeed, 0(/i = — oo) = 1 follows from the fact that the amplitude appertaining to the lower sign between 
the radicals in Eq. (|24p vanishes in the limit of /i — > — oo. 

In accordance with general properties of the subcritical SSB bifurcation [To| , the symmetric solution is expected to 
be stable below the bifurcation point [at N < N^^n, see Eq. d^ ]. and unstable above it. The asymmetric branches 
emerging at iV = iVbif are unstable as long as they go backward. In the present case (a = 0), this means they are 
always unstable, as the respective branches in Figs. [31[b,c) never turn forward. All these expectations are completely 
borne out by the stability analysis performed by means of both direct simulations and computation of stability 
eigenvalues, at finite but small values of a (technical details of the procedure are described in the next section). In 
particular, at > Abif the unstable symmetric state spontaneously transforms into a strongly asymmetric breather 
which features irregular oscillations, but remains robust as a whole (quite similar to an example displayed below 
in Fig. [5] for a — 0.7). On the other hand, unstable asymmetric states transform themselves into breathers which 
maintain the original asymmetry of the unstable state, as shown in Fig. 01 

Lastly, all antisymmetric states in the model with delta-functions are unstable too. Their instability is similar to 
that shown below in Fig. El^b) for a = 1, transforming them into strongly asymmetric breathers. As shown in the 
next section, both asymmetric and antisymmetric states may be stabilized in the general model, with finite a. 

C. Concluding remarks concerning the delta-function model 

The existence of the asymmetric states, i.e., the presence of the SSB effect in the present model, can be easily 
explained by the consideration of the above-mentioned limit of /i ^ — oo. Indeed, in this limit, the spatial scale of 
the solution, which is ~ according to Eq. ([13]), is much smaller than the separation between the two delta- 

functions, 2A = 2. Therefore, the full solution effectively splits into a superposition of those independently supported 
by each delta-function in isolation. Further, it is obvious that, for given large Eq. (|12p with an individual 
delta-function gives rise to two solutions: a trivial one, </> = 0, and 



0±(x)=±(2H)^/%xp(-V2R|el), 



(26) 
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FIG. 3: (Color online) Families of asymmetric states in the plane of {O, ji) (a) and {Q, N) (b) for different fixed values of 
nonlinear-potential parameter a [recall that asymmetry ratio Q is defined by Eq. (|25p . solution families with = being 
symmetric ones; panel (b) takes into regard the fact that symmetric states do not exist at A'^ < 1/2]. Plots labeled "a = 0" 
were generated by analytical solutions (|13p . (|14p and ((24}. Panel (c) additionally displays solution branches from (b) near the 
bifurcation point, and (d) shows the coordinates of the bifurcation point, A^'bif and /^bif, as functions of a. Here and in Fig. [TO] 
below, solid and dashed lines depict stable and unstable solutions, respectively. Note that the bifurcation observed in panels 
(b) and (c) at a = is a "fully backward" one: the branches of asymmetric states, which emerge at the bifurcation point, never 
turn forward and, accordingly, always remain unstable. 



where ^ = a;-|-l or^ = 2;— 1; note that the norm of solution is TV = 1, for any fx. The corresponding symmetric 
and antisymmetric states are built, respectively, as superpositions of solutions 0+ (or, equivalently, </>_) centered at 
X = —1 and X = +1, or 0_ centered at x = —1 and (j)^ centered at a; = +1. Asymmetric solutions are represented, in 
the same limit, by a superposition of solution (j)± centered at x — —1 and zero solution around a; = +1, or vice versa. 
Of course, finding the bifurcation point requires one to perform the analysis of the model at finite /i, as done in the 
analytical form above for the case of the delta-functions, and will be done in a numerical form below for the general 
case of finite a in Eq. ([6]). 

It is relevant to compare the above exact results with those which can be easily obtained in the linear counterpart 
of the model, i.e., the one with the DWP based on the set of two delta- functions; as mentioned above, such a linear 
potential may be a plausible ingredient of a more general model, relevant to the description of NDWP settings in 
optics. The stationary version of the linear equation reduces to 



Ai0 + (1/2),?;." -I- e [Six + 1) + S{x -l)]4>^ 0, 



(27) 
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FIG. 4: (Color online) The evolution of an unstable asymmetric state for a = 1/60 and /i — —0.1. 



with constant e > 0. Symmetric solutions, which must be continuous at = 1, are sought for as 

exp (- V2¥l (kl - 1)) , at kl > 1, 
(f>{x) = ' - ^ - - ^ 



sech f yj2 j cosh f yj2 \^i\x \ , at |a;| < 1, 



(28) 



of. Eq. The jump condition, A (0') \x=±i = —2e4>{x — ±1), yields equation •\/2"i7t[ 1 + tanh ^y^Tj/If^ = 2e, 

which has a single solution for at any positive e, i.e., the linear model always supports exactly one symmetric state. 
Antisymmetric solutions are sought for as 



{x) = 



sgn(x) • exp y^^2\ir\ {\x\ - 
(y2p) sinh (^2" 



cosech I 



1)) , at |x| > 1, 
at |x| 



fi\X 



< 1, 



and the respective jump condition takes the form of y/2 1 + coth yy/2 j = 2e. The latter equation has no 

solutions for e < 1/2, and exactly one solution for e > 1/2. Thus, symmetric solution ([28|) is the single state in 
the linear DWP model at e < 1/2, while at e > 1/2 the linear model supports precisely two states, symmetric and 
antisymmetric ones. 

Lastly, the combined model, which includes both the linear potential and its nonlinear counterpart (that may be 
self-attractive, as above, or self-repulsive too, in this case), is also solvable in the case when these features are based 
on the pair of delta-functions. The combined model is described by the following stationary equation, cf. Eqs. 
and dSZl): 

/x0 + (1/2)0" -I- [5{x + 1) + 5{x - 1)] (e0 + (70^) = 0, (29) 

where ct = +1 and —1 corresponds to the nonlinear attraction and repulsion, respectively. In particular, the bifurcation 
occurs only on the branch of the symmetric solutions in the case of tr = -1-1, and only on the antisymmetric branch - 
in the opposite case (self- repulsion). In either case, the value of the chemical potential (/i < 0) at the symmetry- or 
antisymmetry-breaking bifurcation is determined by the following transcendental equation: 



l-e-V2|H 



(30) 



[for e = and a = +1, it reduces to Eq. (HH)]. With a = -1-1, Eq. ([50)1 has exactly one solution for any e > —1/4 and 
no solutions for e < — 1/4 (negative e corresponds to competition between the repulsive linear potential and attractive 
nonlinear pseudopotential). With a — —1, Eq. (l30|) has exactly one solution for any e > 3/4, and no solutions for 
e < 3/4. Analysis of the stability of states found in the combined model is beyond the scope of this work, and in the 
rest of the paper we consider the model without the linear potential. 
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IV. NUMERICAL RESULTS FOR THE GENERAL MODEL 



A. Numerical methods 



To construct localized solutions to stationary equation (fTO|) of the symmetric, asymmetric and antisymmetric types, 
the Newton iterative method was used, starting with the following inputs: 

iM=^))sy„m = B sech (2 (x + 1))+A sech (2 {x - 1)) , 



('^o(a;))antisym = ^ ^cch (x) sin (x) , 

with B — A for symmetric solutions. Then, as mentioned above, numerical analysis of the stability of the stationary 
solutions was performed in two different ways: first, by means of direct simulations of the evolution of slightly 
perturbed solutions, and then through computation of (in)stability eigenvalues for modes of small perturbations. In 
the former case, the stability was tested by adding arbitrary perturbations to the initial state, at the level of ^ 1% of 
the amplitude (in particular, care was taken to test effects of perturbations whose symmetry is different from that of 
the stationary state, such as antisymmetric perturbations added to symmetric states, and vice versa). 
For the computation of eigenvalues, perturbed solutions were looked for as 



iPix, t) = e-'^* |0(a;) + 77 [m (x) e*^* + v* (x) e*^**] } , (31) 



where (j}{x) is a stationary solution to Eq. (|10[) with chemical potential /i, while u and v are components of a 
perturbation mode with an infinitesimal amplitude 77, pertaining to instability growth rate A = Ar + iXi. The 
substitution of expression ((3T|) into Eq. ([5]) and linearization lead to the eigenvalue problem based on the following 
equations: 



V 



2 dx^ 

-g{x)4>\x) +i^+^_2.g(x)02(:,) 




(32) 



The underlying solution, 4>{x), is stable if all eigenvalues associated with it have Ai = 0. Equations ([5^ were solved 
numerically with the help of a finite-difference method. Conclusions concerning the stability of the patterns, drawn 
from direct simulations, always complied with results produced by the computation of eigenvalue. 



B. Results 



The first significant change against the results reported above for the model with the delta- functions (a = 0), which 
happens with the increase of a, is quick stabilization of asymmetric states with larger values of the norm, while ones 
with smaller N remain unstable, originally. At a > 0.2, the symmetric states are stable for all values of N at which 
they exist. Another notable feature of the bifurcation diagrams at finite a, demonstrated by Figs. [2Ib) and[3](d), is 
that the norm at which the SSB bifurcation takes place, iVbif, first decreases with the growth of a from small values 
up to a « 0.8, and then increases with the further growth of a. 

Close to the their stabilization threshold (in particular, at a = 0.2), asymmetric states with a smaller norm, which 
are still unstable, demonstrate a scenario of the instability development different from what was shown by their 
counterpart in Fig. [4] in the case of very small a. Namely, slow regular oscillations, observed in Fig. [5] in this case, 
imply a dynamical re-symmetrization of the imstable asymmetric state. Indeed, densities ["(/"(a;)!^, taken at points 
a; = ±1, perform identical periodic oscillations, with a phase shift of tt between them, as shown in[SJb). 

The stabilization of the asymmetric states at small finite values of a is explained by the change in the character of the 
SSB bifurcation: at a 7^ 0, there appear turning points on branches of asymmetric solutions in the bifurcation diagram, 
cf. Fig. ^c). Past the turning point, the branch goes forward as a stable one. In fact. Fig. ^c) demonstrates a quick 
transformation, with the increase of a, of the subcritical bifurcation into a supercritical one. When the bifurcation is 
supercritical, branches of the asymmetric solutions emerge as stable ones at the bifurcation point, and immediately 
go forward. 

We do not display the quick transition from the sub- to supercritical bifurcation in full detail, as it actually happens 
at very small a, in the range of a < 0.1. The physical estimates given in Section II suggest that so small values of 
the scaled width of the nonlinear-potential wells correspond to physical widths < 1 fim. It seems doubtful that the 
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FIG. 5: (Color online) (a) The evolution of a weakly unstable asymmetric state with a — 0.2 for /i — —0.075 (very close to the 
bifurcation point and stabilization threshold), (b) For the same case, solid and dotted curves show densities j-)/' {x = —1)1'^ and 
[ip (x ~ +1)1^ as functions of time. 



Feshbach-resonance technique would allow one to create a strong local inhomogeneity of the scattering length on 
such a small scale (nevertheless, the exact analytical solutions obtained for a = 0, which provide clear clues for the 
understanding of the general model, are definitely relevant). An additional problem impeding the full analysis of 
the case of very small a is that, in this case, the accumulation of systematic numerical results requires very heavy 
simulations, as the stepsize of the spatial grid must be made much smaller than a. 

Above the bifurcation point, symmetric states found at finite a demonstrate the familiar SSB instability, sponta- 
neously transforming themselves into slightly nonstationary robust modes (breathers), quite close in their shape to 
respective stable asymmetric solitons. A typical example of this transformation in displayed in Fig. [6l 




FIG. 6: The evolution of an unstable symmetric state, at a = 0.7, n — —2.677 and = 10. 



As concerns antisymmetric solutions, both stable and unstable ones have been found at finite a, as illustrated by 
Figs. [7] and [HI Panel (b) of the former figure shows that the density profile of unstable antisymmetric states evolves 
from the double-peak pattern into an asymmetric single-peak one, which features persistent intrinsic oscillations. 
This outcome of the instability development complies with the fact that the instability of the antisymmetric states is 
oscillatory, being accounted for by a quartet of eigenvalues, as seen in Fig. lU^b). In other words, the transition from 
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stable to unstable antisymmetric states may be considered as the Hamiltonian Hopf bifurcation [41 




(a) (b) 

FIG. 7: (Color online) (a) An example of a stable antisymmetric state with a = 1, for = —1.3 and = 6.7137. (b) The 
evolution of an unstable antisymmetric state, also with a = 1, but for = —1.5 and iV = 7.1273. 




X X x X X X 

r r r 

(a) (b) (c) 

FIG. 8: Examples of stable (a,c) and unstable (b) antisymmetric states, found with a fixed norm, A'^ = 10. In each subplot, 
the left and right panels show, severally, profiles of the stationary states and spectral planes of the (in)stability eigenvalues. 
Parameters are (a) a = 0.80 and /i — —3.05; (b) a — 1.11 and /i — —3.14; (c) a — 1.30 and fi — —3.07. 

Figure [9] displays a combined diagram in the plane of the norm of the solution and width of the nonlinear potential 
wells, N and a, which summarizes the existence and stability results for the states of all the three types - symmetric, 
asymmetric and antisymmetric ones. The dashed-dotted line in the figure designates the symmetry-breaking bifurca- 
tion. Solely symmetric states exist below this line (they are stable in that region), and stable asymmetric states exist 
above the line, where the symmetric ones are unstable. Solid curves in Fig. [9]depict stability borders of antisymmetric 
solutions. 

For the reasons explained above, the region of very small values of a, where the "quick" stabilization of asymmetric 
states takes place, is not included. However, the region of the existence of the analytical symmetric and antisymmetric 
solutions in the model with delta-functions (a = 0), and the respective bifurcation point, as given by Eq. (|23|) . are 
shown by the bold vertical segment and square-marked dot on the axis of a = (recall that the exact asymmetric 
solutions are unstable above the bifurcation point in the model with a = 0). 

Because the modulation profile jl]) does not feature the double- well structure for a > \/2 (see Fig. [TJ, the SSB 
bifurcation tends to disappear as a approaches \/2. Actual results are included in Fig. [9] for a < 1.35, as the 
convergence of the numerical scheme becomes poor for values of a still closer to v^- 

The chain of circles in Figure [5] designates the threshold (minimum norm), A^min, necessary for the existence of 
symmetric states in the model. As mentioned above, in the case of a = the exact threshold is Ny^inia = 0) = 1/2, 
and it is observed in Fig. [5] that the threshold remains in the ballpark of this value at finite a, which can be easily 
explained. Indeed, the minimum of N is attained at /i ^ —0, in which limit the spatial scale of the wave function, 
~ l/-\/2|/i|, is much larger than the size of the NDWP structure, 2A = 2. Thus, from the viewpoint of this weakly 
localized wave function, modulation pattern ([6]) looks like 2S(x). In the corresponding approximation, the wave 
function takes the form of expression ((26)) divided by \/2, and the respective norm is, indeed, 1/2. 

For TV < A'min, the condensate confined to the trap of large length L (i.e., in the thermodynamic limit) will tend 
to form a quasi- uniform nearly linear state, with 4>{x) = ^/N/L. As follows from Eqs. ([5]) and the energy of the 
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FIG. 9: (Color online) Stability and existence borders for symmetric, asymmetric and antisymmetric states in the plane of 
{a,N). The vertical dashed-dotted line corresponds to a = \/2, beyond which the double- well structure in Eq. ^ turns into 
a single-well one. The lower dashed-dotted line depicts the symmetry-breaking bifurcation, with stable asymmetric solutions 
existing above it. The chain of circles in the bottom of the parameter plane designates the minimum norm necessary for the 
existence of symmetric states, which are stable between the existence border and the bifurcation line. Label "Stable Asym - 
Anti" marks regions of the bistability, where both asymmetric and antisymmetric states are stable. 



small-amplitude uniform state is 

Ho ~ -N^/L\ (33) 

In fact, this state realizes a minimum of the energy (cf. Fig. \TU\i . i.e., the system's ground state. Nevertheless, a well- 
known fact is that dynamically stable localized states different from the ground state, such as the above-mentioned 
gap solitons in the repulsive condensate [30j . or their broader counterparts, in the form of the so-called gap- waves 
|49|, can be created in the experiment. 

Another notable feature observed in Fig. [9] is the bistability, i.e., the coexistence of stable asymmetric and an- 
tisymmetric states above the stability border of the latter state. In fact, the bistability always takes place when 
antisymmetric states are stable. It is interesting too that the stability area for the antisymmetric states consists of 
two separate regions. Finally, it is relevant to mention that, as well as in the analytically solvable model with the 
delta- functions {a —^ 0), the antisymmetric states never undergo a bifurcation at finite a. 

We stress that the stability borders displayed in Fig. [5] were identified by means of direct simulations and the 
computation of stability eigenvalues, both methods yielding identical results. In particular, the sets of eigenvalues 
displayed in Fig. [5] clearly confirm the presence of two disjoint stability areas for antisymmetric states. 

In the case of the bistability involving the asymmetric and antisymmetric states, it is interesting to compare their 
energies (values of the Hamiltonian) . To this end. Fig. [TO] displays a typical example of the dependence of H on norm 
iV. The situation observed in this figure is also true in the general case: stable antisymmetric states realize smaller 
values of H than the asymmetric counterparts coexisting with them. However, as argued above, dynamically stable 
states can be created in the experiment even if their energy is higher than in some competing states. In particular, 
Fig. [S] demonstrates that an unstable symmetric state definitely self-traps into an asymmetric robust breather (which 
is close to a stable stationary solution), despite the fact that a stable antisymmetric state exists at the same values 
of A'^ = 10 and a = 0.7, as seen from Fig. [SI 

Note that all curves in Fig. [TO] start from finite threshold values of N corresponding, as said above, to the minimum 
norm (iVmin) necessary for the existence of the respective states. In particular, for the branch of symmetric solutions, 
.^min is close to 1/2, as argued above (cf. the existence border in the bottom of Fig. while the asymmetric branch 
originates at the bifurcation point (in agreement with the location of the respective dashed-dotted line in Fig. [5]), 
at which the symmetric solution loses its stability. The branch of antisymmetric solutions features a fold in Fig. [TO] 
(in the region where these solutions are unstable), which is similar to the the above-mentioned fact that dependence 
N{fi) in Eq. ([^T|l for the unstable exact antisymmetric states has a minimum, (A^antisym)^^;^ ~ 1-84 at /i w —0.58. If 
replotted in terms of H and N, Eq. ([^ features a similar fold, at iV = (A^antisym)jjjiij- 
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FIG. 10: (color online) The Hamiltonian versus the norm for solutions of different types, with fixed a = 0.75. 

V. CONCLUSION 

In this work, we have introduced a model of the nonhnear double- well potential (NDWP), alias a double- well 
pseudopotential, which can be created in BEC, by means of the spatially inhomogeneous Feshbach resonance, and 
also in nonlinear optics. The model provides for a previoulsy unexplored setting in which effects of the spontaneous 
symmetry breaking (SSB) can be studied. 

In the limit when each potential well is induced by the delta-function, full analytical solutions were obtained for 
symmetric, antisymmetric and asymmetric states. The symmetric states are stable in that case up to the symmetry- 
breaking bifurcation point, but beyond the bifurcation both symmetric and emergent asymmetric states are unstable. 
In particular, the asymmetric configurations transform themselves into breathers. The instability of all the stationary 
asymmetric states in the model with the delta-functions is explained by fact that the respective SSB bifurcation is of 
a "fully backward" type, with branches of the asymmetric solutions never turning forward. All antisymmetric states 
are unstable too, in this limit form of the model. 

The increase of the width of the potential wells readily stabilizes the asymmetric states, which concurs with the 
change of the character of the SSB bifurcation from sub- to supercritical. Close to the stabilization border, unstable 
asymmetric states develop slow intrinsic oscillations, featuring effective dynamical re-symmetrization. Antisymmetric 
states may also be stable in the NDWP structure with a finite width of the wells, which implies the bistability 
between asymmetric and antisymmetric states. The symmetric states exist above a finite threshold, in terms of the 
norm (number of atoms in the condensate), and they develop the usual SSB instability above the bifurcation point. 
A simple explanation to the existence threshold was given, and an integrated diagram for the existence and stability 
of the trapped states of all the three types has been produced. 

The analysis presented in this work suggests new experiments in the matter-wave and nonlinear-optical settings. 
The analysis can also be developed in other directions. In particular, it may be interesting to study a two-dimensional 
nonlinear-DWP configuration. In the 2D space, a triangular configuration with three nonlinear (pseudo-)potential 
wells may be considered too. 
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the Israel Science Foundation through the Center-of- Excellence grant No. 8006/03, and by the Thailand Research 
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